##Copy the part starting with source() into the console and run it there - this will create a log file that includes the output

#source("~/work/November_2024/November_2024_PanelMatchScript_Robustness_25SHORT.R", echo = T, max.deparse.length=10000)
rm(list=ls(all=TRUE))
sink("~/work/November_2024/R_out/log_files/full_log_mainROBUST.log", split = T)

# ======================================================================
# Project:    Closed borders, closed minds? COVID-related border closures,
#             EU support and hostility towards immigrants
#
# Script:     Appendix Figures D5-D11
#
# Authors:    Lisa Herbig (l.j.herbig@uva.nl)
#             Asli Unan (a.unan@uva.nl)
#
# Date:       2025-03-24
# ======================================================================

# Description:
# This script generates Appendix Figures D5-D11, performing various robustness checks.

# ----------------------------------------------------------------------
# Load libraries
# ----------------------------------------------------------------------
library("did2s")
library("dplyr")
library("foreign")
library("ggplot2")
library("haven")
library("lmtest")
library("lubridate")
library("modelsummary")
library("PanelMatch")
library("panelView")
library("readr")    
library("tidyr")
library("tidyverse")


## Read the files in
df_iso_work_dd = readRDS("~/work/November_2024/R_out/df_iso_work_dd_25.rds")

# Figure D5 and D6 Lockdown control --------------------------------------------------------

#############################################Refugee - Appendix Figure D5
# No Matching
any_closure_binary_control0_nm = PanelMatch(lag = 4, time.id = "weekyear_i",
                                            unit.id = "iso2.y_i", 
                                            treatment = "any_closure_binary_control0", 
                                            refinement.method = "none", 
                                            data = as.data.frame(df_iso_work_dd), qoi = "att" ,
                                            covs.formula = ~ lockdown +East + age_n + female_n + unemployed_n + working_abroad_n + university_n + EU_residence_n  + kr_foreigner_n + CONTACTABROAD_n + HOMEOFFICE_n + Voting_2018_LeftRight_n+ mode_phone_n +  iso2_pop_n + I(lag(kr_inz_rate_n, 1:4)),
                                            outcome.var = "immigration_pca_n",
                                            size.match = 1,
                                            match.missing = FALSE,
                                            matching = FALSE,
                                            lead = 1:4, forbid.treatment.reversal = FALSE, 
                                            placebo.test = T)

pt_DL_NN <- placebo_test(any_closure_binary_control0_nm, data = as.data.frame(df_iso_work_dd), plot = F, se.method = "conditional")
pt_DL_NN_t_4_est <- pt_DL_NN$estimates[1]
pt_DL_NN_t_3_est <- pt_DL_NN$estimates[2]
pt_DL_NN_t_2_est <- pt_DL_NN$estimates[3]

pt_DL_NN_t_4_high <- pt_DL_NN_t_4_est+(pt_DL_NN$standard.errors[1]*1.96)
pt_DL_NN_t_3_high <- pt_DL_NN_t_3_est+(pt_DL_NN$standard.errors[2]*1.96)
pt_DL_NN_t_2_high <- pt_DL_NN_t_2_est+(pt_DL_NN$standard.errors[3]*1.96)
pt_DL_NN_t_4_low <- pt_DL_NN_t_4_est-(pt_DL_NN$standard.errors[1]*1.96)
pt_DL_NN_t_3_low <- pt_DL_NN_t_3_est-(pt_DL_NN$standard.errors[2]*1.96)
pt_DL_NN_t_2_low <- pt_DL_NN_t_2_est-(pt_DL_NN$standard.errors[3]*1.96)

Results_DL_NN <- PanelEstimate(sets = any_closure_binary_control0_nm, data = as.data.frame(df_iso_work_dd), se.method = "conditional")
summary(Results_DL_NN)

# Matching
any_closure_binary_control0_m = PanelMatch(lag = 4, time.id = "weekyear_i",
                                           unit.id = "iso2.y_i", 
                                           treatment = "any_closure_binary_control0", 
                                           refinement.method = "mahalanobis", 
                                           data = as.data.frame(df_iso_work_dd), qoi = "att",
                                           covs.formula = ~ lockdown +East + age_n + female_n + unemployed_n + working_abroad_n + university_n + EU_residence_n  + kr_foreigner_n + CONTACTABROAD_n + HOMEOFFICE_n + Voting_2018_LeftRight_n + mode_phone_n +  iso2_pop_n + I(lag(kr_inz_rate_n, 1:4)),
                                           outcome.var = "immigration_pca_n",
                                           size.match = 1,
                                           match.missing = TRUE,
                                           matching = TRUE,
                                           lead = 1:4, forbid.treatment.reversal = FALSE, 
                                           placebo.test = T)


pt_DL_mahalanobis <- placebo_test(any_closure_binary_control0_m, data = as.data.frame(df_iso_work_dd), plot = F)

pt_DL_mahalanobis_t_4_est <- pt_DL_mahalanobis$estimates[1]
pt_DL_mahalanobis_t_3_est <- pt_DL_mahalanobis$estimates[2]
pt_DL_mahalanobis_t_2_est <- pt_DL_mahalanobis$estimates[3]

pt_DL_mahalanobis_t_4_high <- pt_DL_mahalanobis_t_4_est+(pt_DL_mahalanobis$standard.errors[1]*1.96)
pt_DL_mahalanobis_t_3_high <- pt_DL_mahalanobis_t_3_est+(pt_DL_mahalanobis$standard.errors[2]*1.96)
pt_DL_mahalanobis_t_2_high <- pt_DL_mahalanobis_t_2_est+(pt_DL_mahalanobis$standard.errors[3]*1.96)
pt_DL_mahalanobis_t_4_low <- pt_DL_mahalanobis_t_4_est-(pt_DL_mahalanobis$standard.errors[1]*1.96)
pt_DL_mahalanobis_t_3_low <- pt_DL_mahalanobis_t_3_est-(pt_DL_mahalanobis$standard.errors[2]*1.96)
pt_DL_mahalanobis_t_2_low <- pt_DL_mahalanobis_t_2_est-(pt_DL_mahalanobis$standard.errors[3]*1.96)

Results_DL_mahalanobis <- PanelEstimate(sets = any_closure_binary_control0_m, data = as.data.frame(df_iso_work_dd), se.method = "conditional")
summary(Results_DL_mahalanobis)

# Generate Figure 2

plot_table_NN <- rbind(pt_DL_NN_t_4_est, pt_DL_NN_t_3_est, pt_DL_NN_t_2_est, 0, as.data.frame(Results_DL_NN$estimates))
colnames(plot_table_NN) <- "estimate"
plot_table_NN$high <- c(pt_DL_NN_t_4_high, pt_DL_NN_t_3_high, pt_DL_NN_t_2_high, 0, Results_DL_NN$estimates + (Results_DL_NN$standard.error*1.96))
plot_table_NN$low <- c(pt_DL_NN_t_4_low, pt_DL_NN_t_3_low, pt_DL_NN_t_2_low, 0, Results_DL_NN$estimates - (Results_DL_NN$standard.error*1.96))
plot_table_NN$Matching <- "No Matching" 
rownames(plot_table_NN)[1:5] <- c("t-4","t-3","t-2","t-1","t")
plot_table_NN$Week <- rownames(plot_table_NN)

plot_table_R <- rbind(pt_DL_mahalanobis_t_4_est, pt_DL_mahalanobis_t_3_est, pt_DL_mahalanobis_t_2_est, 0, as.data.frame(Results_DL_mahalanobis$estimates))
colnames(plot_table_R) <- "estimate"
plot_table_R$high <- c(pt_DL_mahalanobis_t_4_high, pt_DL_mahalanobis_t_3_high, pt_DL_mahalanobis_t_2_high, 0, Results_DL_mahalanobis$estimates + (Results_DL_mahalanobis$standard.error*1.96))
plot_table_R$low <- c(pt_DL_mahalanobis_t_4_low, pt_DL_mahalanobis_t_3_low, pt_DL_mahalanobis_t_2_low, 0, Results_DL_mahalanobis$estimates - (Results_DL_mahalanobis$standard.error*1.96))
plot_table_R$Matching <- "Matching + Refinement" 
rownames(plot_table_R)[1:5] <- c("t-4","t-3","t-2","t-1","t")
plot_table_R$Week <- rownames(plot_table_R)

plot_table_DL <- rbind(plot_table_R, plot_table_NN)
plot_table_DL$Week <- factor(plot_table_DL$Week, levels = c("t-4", "t-3", "t-2", "t-1", "t", "t+1", "t+2", "t+3"))

plot_table_DL$`Estimation Strategy` <- factor(plot_table_DL$Matching, levels = c("No Matching", "Matching + Refinement"))

# Plot
pd = position_dodge(.6)   
ggplot(plot_table_DL, aes(x = Week, y = estimate, color = `Estimation Strategy`)) + geom_hline(yintercept = 0, color = "gray40", linetype="dashed") +
  geom_point(shape = 16, size  = 6, position = pd) +
  geom_errorbar(aes(ymin  = low, ymax  = high), width = 0, size  = 1.45, position = pd) +
  theme_bw() + theme() + ylab("Estimated Treatment Effect on Attitudes towards Refugees") + scale_color_manual(values= c("gray70", "#377EB8")) + theme(legend.text=element_text(size=20), axis.text = element_text(size = 20), axis.title = element_text(size = 20), text = element_text(size = 20), plot.title = element_text(hjust = 0.5)) + ggtitle("Border Closures") 

ggsave(path = "~/work/November_2024/R_out/robustness_plots/", filename = "ROBUST_Lockdown_Refugees_MainPlot.pdf")


#############################################CONNECTION EUROPE - Figure D6
# No Matching
any_closure_binary_control0_nm = PanelMatch(lag = 4, time.id = "weekyear_i",
                                            unit.id = "iso2.y_i", 
                                            treatment = "any_closure_binary_control0", 
                                            refinement.method = "none", 
                                            data = as.data.frame(df_iso_work_dd), qoi = "att" ,
                                            covs.formula = ~ lockdown + East + age_n + female_n + unemployed_n + working_abroad_n + university_n + EU_residence_n  + kr_foreigner_n + CONTACTABROAD_n + HOMEOFFICE_n + Voting_2018_LeftRight_n+ mode_phone_n +  iso2_pop_n + I(lag(kr_inz_rate_n, 1:4)),
                                            outcome.var = "con_europe_n",
                                            size.match = 1,
                                            match.missing = FALSE,
                                            matching = FALSE,
                                            lead = 1:4, forbid.treatment.reversal = FALSE, 
                                            placebo.test = T)

pt_DL_NN <- placebo_test(any_closure_binary_control0_nm, data = as.data.frame(df_iso_work_dd), plot = F, se.method = "conditional")
pt_DL_NN_t_4_est <- pt_DL_NN$estimates[1]
pt_DL_NN_t_3_est <- pt_DL_NN$estimates[2]
pt_DL_NN_t_2_est <- pt_DL_NN$estimates[3]

pt_DL_NN_t_4_high <- pt_DL_NN_t_4_est+(pt_DL_NN$standard.errors[1]*1.96)
pt_DL_NN_t_3_high <- pt_DL_NN_t_3_est+(pt_DL_NN$standard.errors[2]*1.96)
pt_DL_NN_t_2_high <- pt_DL_NN_t_2_est+(pt_DL_NN$standard.errors[3]*1.96)
pt_DL_NN_t_4_low <- pt_DL_NN_t_4_est-(pt_DL_NN$standard.errors[1]*1.96)
pt_DL_NN_t_3_low <- pt_DL_NN_t_3_est-(pt_DL_NN$standard.errors[2]*1.96)
pt_DL_NN_t_2_low <- pt_DL_NN_t_2_est-(pt_DL_NN$standard.errors[3]*1.96)

Results_DL_NN <- PanelEstimate(sets = any_closure_binary_control0_nm, data = as.data.frame(df_iso_work_dd), se.method = "conditional")
summary(Results_DL_NN)

# Matching
any_closure_binary_control0_m = PanelMatch(lag = 4, time.id = "weekyear_i",
                                           unit.id = "iso2.y_i", 
                                           treatment = "any_closure_binary_control0", 
                                           refinement.method = "mahalanobis", 
                                           data = as.data.frame(df_iso_work_dd), qoi = "att",
                                           covs.formula = ~ lockdown +East + age_n + female_n + unemployed_n + working_abroad_n + university_n + EU_residence_n  + kr_foreigner_n + CONTACTABROAD_n + HOMEOFFICE_n + Voting_2018_LeftRight_n + mode_phone_n +  iso2_pop_n + I(lag(kr_inz_rate_n, 1:4)),
                                           outcome.var = "con_europe_n",
                                           size.match = 1,
                                           match.missing = TRUE,
                                           matching = TRUE,
                                           lead = 1:4, forbid.treatment.reversal = FALSE, 
                                           placebo.test = T)


pt_DL_mahalanobis <- placebo_test(any_closure_binary_control0_m, data = as.data.frame(df_iso_work_dd), plot = F)

pt_DL_mahalanobis_t_4_est <- pt_DL_mahalanobis$estimates[1]
pt_DL_mahalanobis_t_3_est <- pt_DL_mahalanobis$estimates[2]
pt_DL_mahalanobis_t_2_est <- pt_DL_mahalanobis$estimates[3]



pt_DL_mahalanobis_t_4_high <- pt_DL_mahalanobis_t_4_est+(pt_DL_mahalanobis$standard.errors[1]*1.96)
pt_DL_mahalanobis_t_3_high <- pt_DL_mahalanobis_t_3_est+(pt_DL_mahalanobis$standard.errors[2]*1.96)
pt_DL_mahalanobis_t_2_high <- pt_DL_mahalanobis_t_2_est+(pt_DL_mahalanobis$standard.errors[3]*1.96)
pt_DL_mahalanobis_t_4_low <- pt_DL_mahalanobis_t_4_est-(pt_DL_mahalanobis$standard.errors[1]*1.96)
pt_DL_mahalanobis_t_3_low <- pt_DL_mahalanobis_t_3_est-(pt_DL_mahalanobis$standard.errors[2]*1.96)
pt_DL_mahalanobis_t_2_low <- pt_DL_mahalanobis_t_2_est-(pt_DL_mahalanobis$standard.errors[3]*1.96)

Results_DL_mahalanobis <- PanelEstimate(sets = any_closure_binary_control0_m, data = as.data.frame(df_iso_work_dd), se.method = "conditional")
summary(Results_DL_mahalanobis)

# Generate Figure 2

plot_table_NN <- rbind(pt_DL_NN_t_4_est, pt_DL_NN_t_3_est, pt_DL_NN_t_2_est, 0, as.data.frame(Results_DL_NN$estimates))
colnames(plot_table_NN) <- "estimate"
plot_table_NN$high <- c(pt_DL_NN_t_4_high, pt_DL_NN_t_3_high, pt_DL_NN_t_2_high, 0, Results_DL_NN$estimates + (Results_DL_NN$standard.error*1.96))
plot_table_NN$low <- c(pt_DL_NN_t_4_low, pt_DL_NN_t_3_low, pt_DL_NN_t_2_low, 0, Results_DL_NN$estimates - (Results_DL_NN$standard.error*1.96))
plot_table_NN$Matching <- "No Matching" 
rownames(plot_table_NN)[1:5] <- c("t-4","t-3","t-2","t-1","t")
plot_table_NN$Week <- rownames(plot_table_NN)

plot_table_R <- rbind(pt_DL_mahalanobis_t_4_est, pt_DL_mahalanobis_t_3_est, pt_DL_mahalanobis_t_2_est, 0, as.data.frame(Results_DL_mahalanobis$estimates))
colnames(plot_table_R) <- "estimate"
plot_table_R$high <- c(pt_DL_mahalanobis_t_4_high, pt_DL_mahalanobis_t_3_high, pt_DL_mahalanobis_t_2_high, 0, Results_DL_mahalanobis$estimates + (Results_DL_mahalanobis$standard.error*1.96))
plot_table_R$low <- c(pt_DL_mahalanobis_t_4_low, pt_DL_mahalanobis_t_3_low, pt_DL_mahalanobis_t_2_low, 0, Results_DL_mahalanobis$estimates - (Results_DL_mahalanobis$standard.error*1.96))
plot_table_R$Matching <- "Matching + Refinement" 
rownames(plot_table_R)[1:5] <- c("t-4","t-3","t-2","t-1","t")
plot_table_R$Week <- rownames(plot_table_R)

plot_table_DL <- rbind(plot_table_R, plot_table_NN)
plot_table_DL$Week <- factor(plot_table_DL$Week, levels = c("t-4", "t-3", "t-2", "t-1", "t", "t+1", "t+2", "t+3"))

plot_table_DL$`Estimation Strategy` <- factor(plot_table_DL$Matching, levels = c("No Matching", "Matching + Refinement"))

# Plot
pd = position_dodge(.6)   
ggplot(plot_table_DL, aes(x = Week, y = estimate, color = `Estimation Strategy`)) + geom_hline(yintercept = 0, color = "gray40", linetype="dashed") +
  geom_point(shape = 16, size  = 6, position = pd) +
  geom_errorbar(aes(ymin  = low, ymax  = high), width = 0, size  = 1.45, position = pd) +
  theme_bw() + theme() + ylab("Estimated Treatment Effect on Attachment with Europe") + scale_color_manual(values= c("gray70", "#377EB8")) + theme(legend.text=element_text(size=20), axis.text = element_text(size = 20), axis.title = element_text(size = 20), text = element_text(size = 20), plot.title = element_text(hjust = 0.5)) + ggtitle("Border Closures") 

ggsave(path = "~/work/November_2024/R_out/robustness_plots/", filename = "ROBUST_Lockdown_Con-Europe_MainPlot.pdf")


#############################################Refugee CULTURE
# No Matching
any_closure_binary_control0_nm = PanelMatch(lag = 4, time.id = "weekyear_i",
                                            unit.id = "iso2.y_i", 
                                            treatment = "any_closure_binary_control0", 
                                            refinement.method = "none", 
                                            data = as.data.frame(df_iso_work_dd), qoi = "att" ,
                                            covs.formula = ~ lockdown +East + age_n + female_n + unemployed_n + working_abroad_n + university_n + EU_residence_n  + kr_foreigner_n + CONTACTABROAD_n + HOMEOFFICE_n + Voting_2018_LeftRight_n+ mode_phone_n +  iso2_pop_n + I(lag(kr_inz_rate_n, 1:4)),
                                            outcome.var = "refugee_culture_n",
                                            size.match = 1,
                                            match.missing = FALSE,
                                            matching = FALSE,
                                            lead = 1:4, forbid.treatment.reversal = FALSE, 
                                            placebo.test = T)

pt_DL_NN <- placebo_test(any_closure_binary_control0_nm, data = as.data.frame(df_iso_work_dd), plot = F, se.method = "conditional")
pt_DL_NN_t_4_est <- pt_DL_NN$estimates[1]
pt_DL_NN_t_3_est <- pt_DL_NN$estimates[2]
pt_DL_NN_t_2_est <- pt_DL_NN$estimates[3]

pt_DL_NN_t_4_high <- pt_DL_NN_t_4_est+(pt_DL_NN$standard.errors[1]*1.96)
pt_DL_NN_t_3_high <- pt_DL_NN_t_3_est+(pt_DL_NN$standard.errors[2]*1.96)
pt_DL_NN_t_2_high <- pt_DL_NN_t_2_est+(pt_DL_NN$standard.errors[3]*1.96)
pt_DL_NN_t_4_low <- pt_DL_NN_t_4_est-(pt_DL_NN$standard.errors[1]*1.96)
pt_DL_NN_t_3_low <- pt_DL_NN_t_3_est-(pt_DL_NN$standard.errors[2]*1.96)
pt_DL_NN_t_2_low <- pt_DL_NN_t_2_est-(pt_DL_NN$standard.errors[3]*1.96)

Results_DL_NN <- PanelEstimate(sets = any_closure_binary_control0_nm, data = as.data.frame(df_iso_work_dd), se.method = "conditional")
summary(Results_DL_NN)

# Matching
any_closure_binary_control0_m = PanelMatch(lag = 4, time.id = "weekyear_i",
                                           unit.id = "iso2.y_i", 
                                           treatment = "any_closure_binary_control0", 
                                           refinement.method = "mahalanobis", 
                                           data = as.data.frame(df_iso_work_dd), qoi = "att",
                                           covs.formula = ~ lockdown +East + age_n + female_n + unemployed_n + working_abroad_n + university_n + EU_residence_n  + kr_foreigner_n + CONTACTABROAD_n + HOMEOFFICE_n + Voting_2018_LeftRight_n + mode_phone_n +  iso2_pop_n + I(lag(kr_inz_rate_n, 1:4)),
                                           outcome.var = "refugee_culture_n",
                                           size.match = 1,
                                           match.missing = TRUE,
                                           matching = TRUE,
                                           lead = 1:4, forbid.treatment.reversal = FALSE, 
                                           placebo.test = T)


pt_DL_mahalanobis <- placebo_test(any_closure_binary_control0_m, data = as.data.frame(df_iso_work_dd), plot = F)

pt_DL_mahalanobis_t_4_est <- pt_DL_mahalanobis$estimates[1]
pt_DL_mahalanobis_t_3_est <- pt_DL_mahalanobis$estimates[2]
pt_DL_mahalanobis_t_2_est <- pt_DL_mahalanobis$estimates[3]

pt_DL_mahalanobis_t_4_high <- pt_DL_mahalanobis_t_4_est+(pt_DL_mahalanobis$standard.errors[1]*1.96)
pt_DL_mahalanobis_t_3_high <- pt_DL_mahalanobis_t_3_est+(pt_DL_mahalanobis$standard.errors[2]*1.96)
pt_DL_mahalanobis_t_2_high <- pt_DL_mahalanobis_t_2_est+(pt_DL_mahalanobis$standard.errors[3]*1.96)
pt_DL_mahalanobis_t_4_low <- pt_DL_mahalanobis_t_4_est-(pt_DL_mahalanobis$standard.errors[1]*1.96)
pt_DL_mahalanobis_t_3_low <- pt_DL_mahalanobis_t_3_est-(pt_DL_mahalanobis$standard.errors[2]*1.96)
pt_DL_mahalanobis_t_2_low <- pt_DL_mahalanobis_t_2_est-(pt_DL_mahalanobis$standard.errors[3]*1.96)

Results_DL_mahalanobis <- PanelEstimate(sets = any_closure_binary_control0_m, data = as.data.frame(df_iso_work_dd), se.method = "conditional")
summary(Results_DL_mahalanobis)

# Generate Figure 2

plot_table_NN <- rbind(pt_DL_NN_t_4_est, pt_DL_NN_t_3_est, pt_DL_NN_t_2_est, 0, as.data.frame(Results_DL_NN$estimates))
colnames(plot_table_NN) <- "estimate"
plot_table_NN$high <- c(pt_DL_NN_t_4_high, pt_DL_NN_t_3_high, pt_DL_NN_t_2_high, 0, Results_DL_NN$estimates + (Results_DL_NN$standard.error*1.96))
plot_table_NN$low <- c(pt_DL_NN_t_4_low, pt_DL_NN_t_3_low, pt_DL_NN_t_2_low, 0, Results_DL_NN$estimates - (Results_DL_NN$standard.error*1.96))
plot_table_NN$Matching <- "No Matching" 
rownames(plot_table_NN)[1:5] <- c("t-4","t-3","t-2","t-1","t")
plot_table_NN$Week <- rownames(plot_table_NN)

plot_table_R <- rbind(pt_DL_mahalanobis_t_4_est, pt_DL_mahalanobis_t_3_est, pt_DL_mahalanobis_t_2_est, 0, as.data.frame(Results_DL_mahalanobis$estimates))
colnames(plot_table_R) <- "estimate"
plot_table_R$high <- c(pt_DL_mahalanobis_t_4_high, pt_DL_mahalanobis_t_3_high, pt_DL_mahalanobis_t_2_high, 0, Results_DL_mahalanobis$estimates + (Results_DL_mahalanobis$standard.error*1.96))
plot_table_R$low <- c(pt_DL_mahalanobis_t_4_low, pt_DL_mahalanobis_t_3_low, pt_DL_mahalanobis_t_2_low, 0, Results_DL_mahalanobis$estimates - (Results_DL_mahalanobis$standard.error*1.96))
plot_table_R$Matching <- "Matching + Refinement" 
rownames(plot_table_R)[1:5] <- c("t-4","t-3","t-2","t-1","t")
plot_table_R$Week <- rownames(plot_table_R)

plot_table_DL <- rbind(plot_table_R, plot_table_NN)
plot_table_DL$Week <- factor(plot_table_DL$Week, levels = c("t-4", "t-3", "t-2", "t-1", "t", "t+1", "t+2", "t+3", "t+4"))

plot_table_DL$`Estimation Strategy` <- factor(plot_table_DL$Matching, levels = c("No Matching", "Matching + Refinement"))

# Plot
pd = position_dodge(.6)   
ggplot(plot_table_DL, aes(x = Week, y = estimate, color = `Estimation Strategy`)) + geom_hline(yintercept = 0, color = "gray40", linetype="dashed") +
  geom_point(shape = 16, size  = 6, position = pd) +
  geom_errorbar(aes(ymin  = low, ymax  = high), width = 0, size  = 1.45, position = pd) +
  theme_bw() + theme() + ylab("Estimated Treatment Effect on Attitudes towards Refugees (Culture)") + scale_color_manual(values= c("gray70", "#377EB8")) + theme(legend.text=element_text(size=20), axis.text = element_text(size = 20), axis.title = element_text(size = 20), text = element_text(size = 20), plot.title = element_text(hjust = 0.5)) + ggtitle("Border Closures") 

ggsave(path = "~/work/November_2024/R_out/robustness_plots/", filename = "ROBUST_Lockdown_RefugeesCULTURE_MainPlot.pdf")

# Figure D7 and D8 InzRate Neighboring region control ------------------------

#############################################Refugee - Appendix Figure D7
# No Matching
any_closure_binary_control0_nm = PanelMatch(lag = 4, time.id = "weekyear_i",
                                            unit.id = "iso2.y_i", 
                                            treatment = "any_closure_binary_control0", 
                                            refinement.method = "none", 
                                            data = as.data.frame(df_iso_work_dd), qoi = "att" ,
                                            covs.formula = ~ var_rate_14_day_per_100k + East + age_n + female_n + unemployed_n + working_abroad_n + university_n + EU_residence_n  + kr_foreigner_n + CONTACTABROAD_n + HOMEOFFICE_n + Voting_2018_LeftRight_n+ mode_phone_n +  iso2_pop_n + I(lag(kr_inz_rate_n, 1:4)),
                                            outcome.var = "immigration_pca_n",
                                            size.match = 1,
                                            match.missing = FALSE,
                                            matching = FALSE,
                                            lead = 1:4, forbid.treatment.reversal = FALSE, 
                                            placebo.test = T)

pt_DL_NN <- placebo_test(any_closure_binary_control0_nm, data = as.data.frame(df_iso_work_dd), plot = F, se.method = "conditional")
pt_DL_NN_t_4_est <- pt_DL_NN$estimates[1]
pt_DL_NN_t_3_est <- pt_DL_NN$estimates[2]
pt_DL_NN_t_2_est <- pt_DL_NN$estimates[3]

pt_DL_NN_t_4_high <- pt_DL_NN_t_4_est+(pt_DL_NN$standard.errors[1]*1.96)
pt_DL_NN_t_3_high <- pt_DL_NN_t_3_est+(pt_DL_NN$standard.errors[2]*1.96)
pt_DL_NN_t_2_high <- pt_DL_NN_t_2_est+(pt_DL_NN$standard.errors[3]*1.96)
pt_DL_NN_t_4_low <- pt_DL_NN_t_4_est-(pt_DL_NN$standard.errors[1]*1.96)
pt_DL_NN_t_3_low <- pt_DL_NN_t_3_est-(pt_DL_NN$standard.errors[2]*1.96)
pt_DL_NN_t_2_low <- pt_DL_NN_t_2_est-(pt_DL_NN$standard.errors[3]*1.96)

Results_DL_NN <- PanelEstimate(sets = any_closure_binary_control0_nm, data = as.data.frame(df_iso_work_dd), se.method = "conditional")
summary(Results_DL_NN)

# Matching
any_closure_binary_control0_m = PanelMatch(lag = 4, time.id = "weekyear_i",
                                           unit.id = "iso2.y_i", 
                                           treatment = "any_closure_binary_control0", 
                                           refinement.method = "mahalanobis", 
                                           data = as.data.frame(df_iso_work_dd), qoi = "att",
                                           covs.formula = ~ var_rate_14_day_per_100k + East + age_n + female_n + unemployed_n + working_abroad_n + university_n + EU_residence_n  + kr_foreigner_n + CONTACTABROAD_n + HOMEOFFICE_n + Voting_2018_LeftRight_n + mode_phone_n +  iso2_pop_n + I(lag(kr_inz_rate_n, 1:4)),
                                           outcome.var = "immigration_pca_n",
                                           size.match = 1,
                                           match.missing = TRUE,
                                           matching = TRUE,
                                           lead = 1:4, forbid.treatment.reversal = FALSE, 
                                           placebo.test = T)

pt_DL_mahalanobis <- placebo_test(any_closure_binary_control0_m, data = as.data.frame(df_iso_work_dd), plot = F)

pt_DL_mahalanobis_t_4_est <- pt_DL_mahalanobis$estimates[1]
pt_DL_mahalanobis_t_3_est <- pt_DL_mahalanobis$estimates[2]
pt_DL_mahalanobis_t_2_est <- pt_DL_mahalanobis$estimates[3]

pt_DL_mahalanobis_t_4_high <- pt_DL_mahalanobis_t_4_est+(pt_DL_mahalanobis$standard.errors[1]*1.96)
pt_DL_mahalanobis_t_3_high <- pt_DL_mahalanobis_t_3_est+(pt_DL_mahalanobis$standard.errors[2]*1.96)
pt_DL_mahalanobis_t_2_high <- pt_DL_mahalanobis_t_2_est+(pt_DL_mahalanobis$standard.errors[3]*1.96)
pt_DL_mahalanobis_t_4_low <- pt_DL_mahalanobis_t_4_est-(pt_DL_mahalanobis$standard.errors[1]*1.96)
pt_DL_mahalanobis_t_3_low <- pt_DL_mahalanobis_t_3_est-(pt_DL_mahalanobis$standard.errors[2]*1.96)
pt_DL_mahalanobis_t_2_low <- pt_DL_mahalanobis_t_2_est-(pt_DL_mahalanobis$standard.errors[3]*1.96)

Results_DL_mahalanobis <- PanelEstimate(sets = any_closure_binary_control0_m, data = as.data.frame(df_iso_work_dd), se.method = "conditional")
summary(Results_DL_mahalanobis)

# Generate Figure 2

plot_table_NN <- rbind(pt_DL_NN_t_4_est, pt_DL_NN_t_3_est, pt_DL_NN_t_2_est, 0, as.data.frame(Results_DL_NN$estimates))
colnames(plot_table_NN) <- "estimate"
plot_table_NN$high <- c(pt_DL_NN_t_4_high, pt_DL_NN_t_3_high, pt_DL_NN_t_2_high, 0, Results_DL_NN$estimates + (Results_DL_NN$standard.error*1.96))
plot_table_NN$low <- c(pt_DL_NN_t_4_low, pt_DL_NN_t_3_low, pt_DL_NN_t_2_low, 0, Results_DL_NN$estimates - (Results_DL_NN$standard.error*1.96))
plot_table_NN$Matching <- "No Matching" 
rownames(plot_table_NN)[1:5] <- c("t-4","t-3","t-2","t-1","t")
plot_table_NN$Week <- rownames(plot_table_NN)

plot_table_R <- rbind(pt_DL_mahalanobis_t_4_est, pt_DL_mahalanobis_t_3_est, pt_DL_mahalanobis_t_2_est, 0, as.data.frame(Results_DL_mahalanobis$estimates))
colnames(plot_table_R) <- "estimate"
plot_table_R$high <- c(pt_DL_mahalanobis_t_4_high, pt_DL_mahalanobis_t_3_high, pt_DL_mahalanobis_t_2_high, 0, Results_DL_mahalanobis$estimates + (Results_DL_mahalanobis$standard.error*1.96))
plot_table_R$low <- c(pt_DL_mahalanobis_t_4_low, pt_DL_mahalanobis_t_3_low, pt_DL_mahalanobis_t_2_low, 0, Results_DL_mahalanobis$estimates - (Results_DL_mahalanobis$standard.error*1.96))
plot_table_R$Matching <- "Matching + Refinement" 
rownames(plot_table_R)[1:5] <- c("t-4","t-3","t-2","t-1","t")
plot_table_R$Week <- rownames(plot_table_R)

plot_table_DL <- rbind(plot_table_R, plot_table_NN)
plot_table_DL$Week <- factor(plot_table_DL$Week, levels = c("t-4", "t-3", "t-2", "t-1", "t", "t+1", "t+2", "t+3"))

plot_table_DL$`Estimation Strategy` <- factor(plot_table_DL$Matching, levels = c("No Matching", "Matching + Refinement"))

# Plot
pd = position_dodge(.6)   
ggplot(plot_table_DL, aes(x = Week, y = estimate, color = `Estimation Strategy`)) + geom_hline(yintercept = 0, color = "gray40", linetype="dashed") +
  geom_point(shape = 16, size  = 6, position = pd) +
  geom_errorbar(aes(ymin  = low, ymax  = high), width = 0, size  = 1.45, position = pd) +
  theme_bw() + theme() + ylab("Estimated Treatment Effect on Attitudes towards Refugees") + scale_color_manual(values= c("gray70", "#377EB8")) + theme(legend.text=element_text(size=20), axis.text = element_text(size = 20), axis.title = element_text(size = 20), text = element_text(size = 20), plot.title = element_text(hjust = 0.5)) + ggtitle("Border Closures") 

ggsave(path = "~/work/November_2024/R_out/robustness_plots/", filename = "ROBUST_InNeigh_Refugees_MainPlot.pdf")


#############################################CONNECTION EUROPE - Appendix Figure D8
# No Matching
any_closure_binary_control0_nm = PanelMatch(lag = 4, time.id = "weekyear_i",
                                            unit.id = "iso2.y_i", 
                                            treatment = "any_closure_binary_control0", 
                                            refinement.method = "none", 
                                            data = as.data.frame(df_iso_work_dd), qoi = "att" ,
                                            covs.formula = ~ var_rate_14_day_per_100k + East + age_n + female_n + unemployed_n + working_abroad_n + university_n + EU_residence_n  + kr_foreigner_n + CONTACTABROAD_n + HOMEOFFICE_n + Voting_2018_LeftRight_n+ mode_phone_n +  iso2_pop_n + I(lag(kr_inz_rate_n, 1:4)),
                                            outcome.var = "con_europe_n",
                                            size.match = 1,
                                            match.missing = FALSE,
                                            matching = FALSE,
                                            lead = 1:4, forbid.treatment.reversal = FALSE, 
                                            placebo.test = T)

pt_DL_NN <- placebo_test(any_closure_binary_control0_nm, data = as.data.frame(df_iso_work_dd), plot = F, se.method = "conditional")
pt_DL_NN_t_4_est <- pt_DL_NN$estimates[1]
pt_DL_NN_t_3_est <- pt_DL_NN$estimates[2]
pt_DL_NN_t_2_est <- pt_DL_NN$estimates[3]

pt_DL_NN_t_4_high <- pt_DL_NN_t_4_est+(pt_DL_NN$standard.errors[1]*1.96)
pt_DL_NN_t_3_high <- pt_DL_NN_t_3_est+(pt_DL_NN$standard.errors[2]*1.96)
pt_DL_NN_t_2_high <- pt_DL_NN_t_2_est+(pt_DL_NN$standard.errors[3]*1.96)
pt_DL_NN_t_4_low <- pt_DL_NN_t_4_est-(pt_DL_NN$standard.errors[1]*1.96)
pt_DL_NN_t_3_low <- pt_DL_NN_t_3_est-(pt_DL_NN$standard.errors[2]*1.96)
pt_DL_NN_t_2_low <- pt_DL_NN_t_2_est-(pt_DL_NN$standard.errors[3]*1.96)

Results_DL_NN <- PanelEstimate(sets = any_closure_binary_control0_nm, data = as.data.frame(df_iso_work_dd), se.method = "conditional")
summary(Results_DL_NN)

# Matching
any_closure_binary_control0_m = PanelMatch(lag = 4, time.id = "weekyear_i",
                                           unit.id = "iso2.y_i", 
                                           treatment = "any_closure_binary_control0", 
                                           refinement.method = "mahalanobis", 
                                           data = as.data.frame(df_iso_work_dd), qoi = "att",
                                           covs.formula = ~ var_rate_14_day_per_100k + East + age_n + female_n + unemployed_n + working_abroad_n + university_n + EU_residence_n  + kr_foreigner_n + CONTACTABROAD_n + HOMEOFFICE_n + Voting_2018_LeftRight_n + mode_phone_n +  iso2_pop_n + I(lag(kr_inz_rate_n, 1:4)),
                                           outcome.var = "con_europe_n",
                                           size.match = 1,
                                           match.missing = TRUE,
                                           matching = TRUE,
                                           lead = 1:4, forbid.treatment.reversal = FALSE, 
                                           placebo.test = T)


pt_DL_mahalanobis <- placebo_test(any_closure_binary_control0_m, data = as.data.frame(df_iso_work_dd), plot = F)

pt_DL_mahalanobis_t_4_est <- pt_DL_mahalanobis$estimates[1]
pt_DL_mahalanobis_t_3_est <- pt_DL_mahalanobis$estimates[2]
pt_DL_mahalanobis_t_2_est <- pt_DL_mahalanobis$estimates[3]



pt_DL_mahalanobis_t_4_high <- pt_DL_mahalanobis_t_4_est+(pt_DL_mahalanobis$standard.errors[1]*1.96)
pt_DL_mahalanobis_t_3_high <- pt_DL_mahalanobis_t_3_est+(pt_DL_mahalanobis$standard.errors[2]*1.96)
pt_DL_mahalanobis_t_2_high <- pt_DL_mahalanobis_t_2_est+(pt_DL_mahalanobis$standard.errors[3]*1.96)
pt_DL_mahalanobis_t_4_low <- pt_DL_mahalanobis_t_4_est-(pt_DL_mahalanobis$standard.errors[1]*1.96)
pt_DL_mahalanobis_t_3_low <- pt_DL_mahalanobis_t_3_est-(pt_DL_mahalanobis$standard.errors[2]*1.96)
pt_DL_mahalanobis_t_2_low <- pt_DL_mahalanobis_t_2_est-(pt_DL_mahalanobis$standard.errors[3]*1.96)

Results_DL_mahalanobis <- PanelEstimate(sets = any_closure_binary_control0_m, data = as.data.frame(df_iso_work_dd), se.method = "conditional")
summary(Results_DL_mahalanobis)

# Generate Figure 2

plot_table_NN <- rbind(pt_DL_NN_t_4_est, pt_DL_NN_t_3_est, pt_DL_NN_t_2_est, 0, as.data.frame(Results_DL_NN$estimates))
colnames(plot_table_NN) <- "estimate"
plot_table_NN$high <- c(pt_DL_NN_t_4_high, pt_DL_NN_t_3_high, pt_DL_NN_t_2_high, 0, Results_DL_NN$estimates + (Results_DL_NN$standard.error*1.96))
plot_table_NN$low <- c(pt_DL_NN_t_4_low, pt_DL_NN_t_3_low, pt_DL_NN_t_2_low, 0, Results_DL_NN$estimates - (Results_DL_NN$standard.error*1.96))
plot_table_NN$Matching <- "No Matching" 
rownames(plot_table_NN)[1:5] <- c("t-4","t-3","t-2","t-1","t")
plot_table_NN$Week <- rownames(plot_table_NN)

plot_table_R <- rbind(pt_DL_mahalanobis_t_4_est, pt_DL_mahalanobis_t_3_est, pt_DL_mahalanobis_t_2_est, 0, as.data.frame(Results_DL_mahalanobis$estimates))
colnames(plot_table_R) <- "estimate"
plot_table_R$high <- c(pt_DL_mahalanobis_t_4_high, pt_DL_mahalanobis_t_3_high, pt_DL_mahalanobis_t_2_high, 0, Results_DL_mahalanobis$estimates + (Results_DL_mahalanobis$standard.error*1.96))
plot_table_R$low <- c(pt_DL_mahalanobis_t_4_low, pt_DL_mahalanobis_t_3_low, pt_DL_mahalanobis_t_2_low, 0, Results_DL_mahalanobis$estimates - (Results_DL_mahalanobis$standard.error*1.96))
plot_table_R$Matching <- "Matching + Refinement" 
rownames(plot_table_R)[1:5] <- c("t-4","t-3","t-2","t-1","t")
plot_table_R$Week <- rownames(plot_table_R)

plot_table_DL <- rbind(plot_table_R, plot_table_NN)
plot_table_DL$Week <- factor(plot_table_DL$Week, levels = c("t-4", "t-3", "t-2", "t-1", "t", "t+1", "t+2", "t+3"))

plot_table_DL$`Estimation Strategy` <- factor(plot_table_DL$Matching, levels = c("No Matching", "Matching + Refinement"))

# Plot
pd = position_dodge(.6)   
ggplot(plot_table_DL, aes(x = Week, y = estimate, color = `Estimation Strategy`)) + geom_hline(yintercept = 0, color = "gray40", linetype="dashed") +
  geom_point(shape = 16, size  = 6, position = pd) +
  geom_errorbar(aes(ymin  = low, ymax  = high), width = 0, size  = 1.45, position = pd) +
  theme_bw() + theme() + ylab("Estimated Treatment Effect on Attachment with Europe") + scale_color_manual(values= c("gray70", "#377EB8")) + theme(legend.text=element_text(size=20), axis.text = element_text(size = 20), axis.title = element_text(size = 20), text = element_text(size = 20), plot.title = element_text(hjust = 0.5)) + ggtitle("Border Closures") 

ggsave(path = "~/work/November_2024/R_out/robustness_plots/", filename = "ROBUST_InNeigh_Con-Europe_MainPlot.pdf")

# Appendix Figure D9 and D10 No Covariates -----------------------------------------------------------

######Appendix Figure D9
output_folder = "~/work/November_2024/R_out/robustness_plots/"
PM.results <- PanelMatch(lag = 2, time.id = "weekyear_i",
                         unit.id = "iso2.y_i", 
                         treatment = "any_closure_binary_control0", 
                         refinement.method = "none", 
                         data = as.data.frame(df_iso_work_dd), qoi = "att" ,
                         outcome.var = "immigration_pca_n",
                         match.missing = TRUE,
                         matching = TRUE,
                         lead = 1:4, forbid.treatment.reversal = FALSE)
PE.results1 <- PanelEstimate(sets = PM.results, data = as.data.frame(df_iso_work_dd), se.method = "conditional", pooled = FALSE)
pdf(paste0(output_folder, "REBUST_NoCov_IMI_reverse_any_closure_binary_25", ".pdf"))
plot(PE.results1)
summary(PE.results1)
get_covariate_balance(PM.results$att, as.data.frame(df_iso_work_dd), 
                      covariates = c("lockdown" , "var_rate_14_day_per_100k" ,"East" , "age_n" , "female_n" , 
                                     "unemployed_n" , "working_abroad_n" , "university_n" , "EU_residence_n"  , 
                                     "kr_foreigner_n" , "CONTACTABROAD_n" , "HOMEOFFICE_n" , "Voting_2018_LeftRight_n" , 
                                     "mode_phone_n" ,  "iso2_pop_n" , "kr_inz_rate_n"),plot = FALSE)
dev.off()

plot(PE.results1)


######Appendix Figure D10
output_folder = "~/work/November_2024/R_out/robustness_plots/"
PM.results <- PanelMatch(lag = 2, time.id = "weekyear_i",
                         unit.id = "iso2.y_i", 
                         treatment = "any_closure_binary_control0", 
                         refinement.method = "none", 
                         data = as.data.frame(df_iso_work_dd), qoi = "att" ,
                         outcome.var = "con_europe_n",
                         match.missing = TRUE,
                         matching = TRUE,
                         lead = 1:4, forbid.treatment.reversal = FALSE)
PE.results1 <- PanelEstimate(sets = PM.results, data = as.data.frame(df_iso_work_dd), se.method = "conditional", pooled = FALSE)
pdf(paste0(output_folder, "REBUST_NoCov_EU_reverse_any_closure_binary_25", ".pdf"))
plot(PE.results1)
summary(PE.results1)
get_covariate_balance(PM.results$att, as.data.frame(df_iso_work_dd), 
                      covariates = c("lockdown" , "var_rate_14_day_per_100k" ,"East" , "age_n" , "female_n" , 
                                     "unemployed_n" , "working_abroad_n" , "university_n" , "EU_residence_n"  , 
                                     "kr_foreigner_n" , "CONTACTABROAD_n" , "HOMEOFFICE_n" , "Voting_2018_LeftRight_n" , 
                                     "mode_phone_n" ,  "iso2_pop_n" , "kr_inz_rate_n"),plot = FALSE)
dev.off()

plot(PE.results1)


# Appendix Figure D11 Placebo -----------------------------------------------------------

###Placebo - Satisfaction with Democracy
  output_folder = "~/work/November_2024/R_out/robustness_plots/"

PM.results <- PanelMatch(lag = 2, time.id = "weekyear_i",
                                                     unit.id = "iso2.y_i", 
                                                     treatment = "any_closure_binary_control0", 
                                                     refinement.method = "mahalanobis", 
                                                     data = as.data.frame(df_iso_work_dd), qoi = "att" ,
                                                     outcome.var = "satisf_democracy",
                                                     match.missing = TRUE,
                                                     matching = TRUE,
                                                     covs.formula = ~ East + age + HOMEOFFICE + female + unemployed + working_abroad + university + EU_residence  + 
                                                        kr_foreigner + CONTACTABROAD + Voting_2018_LeftRight+ mode_phone + iso2_pop+  I(lag(kr_inz_rate, 1:4)),
                                                     lead = 1:4, forbid.treatment.reversal = FALSE)

PE.results1 <- PanelEstimate(sets = PM.results, data = as.data.frame(df_iso_work_dd), se.method = "conditional", pooled = FALSE)
pdf(paste0(output_folder, "SATDEM_reverse_any_closure_binary_Placebo_ROBUST_25", ".pdf"))

plot(PE.results1)

summary(PE.results1)




sink()